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Abstract 

The numerical simulation of granular systems of even moderate size is a challenging computational problem. 
In most investigations, either Molecular Dynamics or Event-driven Molecular Dynamics is applied. Here 
we show that in certain cases, mainly (but not exclusively) for static granular packings, the Bottom-to-top 
Reconstruction method allows for the efficient simulation of very large systems. We apply the method to heap 
formation, granular flow in a rotating cylinder and to structure formation in nano-powders. We also present 
an efficient implementation of the algorithm in CH — h, including a benchmark. 



Introduction 

"■^The intention of numerical particle simulations of 
granular matter is (at least) two- fold. First, simu- 
Qations are frequently helpful to explain experimental 
i Results that is, to improve our understanding of the 
static and dynamic behavior of granular matter and 
T 7^o improve its corresponding theoretical description. 
i^-Second, yet more important, particle simulations may 
Q^)e helpful for engineers to construct or optimize tech- 
l/Taical devices to handle and process granular materials. 
1 Having a reliable numerical method at hand, the pre- 
^^liction of granular matter behavior is even possible if 
^diere is no theoretical description due to our limited 
^understanding or for more fundamental reasons [11]. 
C^Tjhe main problem in applying particle simulations 
£^is a standard tool in engineering is the enormous 
'c onsumption of computational power (see [32] for an 
j5)verview). The most universal method for particle 
(^simulation is Molecular Dynamics (MD), also called 
CDdscrete Elements Method, where Newton's equation 
^of motion is solved simultaneously for all particles of 
J s nfhe system. The knowledge of the interaction forces 
r^s sufficient, in principle, to simulate any granular as- 
sembly with arbitrary boundary conditions. Unfortu- 
nately, due to the large stiffness of the particle, ex- 
pressed by the Young modulus of the material, MD 
simulations of granular matter require a very small in- 
tegration time step which implies high computational 
costs. At present, for systematic investigations of sys- 
tems over a longer period, the number of particles is 
restricted to significiantly below 10 5 . On the other 
hand, even a relatively small container as used for ex- 
perimental investigations comprises much more parti- 
cles, sometimes billions. 

The simulation may be accelerated considerably by 
assuming that each contact of particles is an instan- 
taneous event, that is, the duration of contacts van- 
ishes and in the entire system at any time only two 
particles are in contact. Under this condition, Event- 



driven MD (EMD) may be applied. Here, instead of 
integrating Newton's equation, the evolution of the 
system is determined by the sequence of particle colli- 
sions, where the postcollision velocities are computed 
as functions of the precollision velocities and the coef- 
ficients of restitution which are functions of the impact 
velocity themselves, see [32]. In the limit of very dilute 
systems of stiff particles, called granular gases, both 
algorithms, MD and EMD, deliver identical results. 
Using desktop computers, at present we can simulate 
about 10 7 particles using EMD, however, the above- 
mentioned condition of EMD does obviously not allow 
for its application as a general simulation tool. EMD, 
in principle, is not suited to simulate systems where 
particles contact each other over a non-vanishing pe- 
riod of time. 

A complementary approach is the Bottom- To- Top Re- 
construction method (BTR). Here the simultaneous 
numerical integration of the iV-particle system is ap- 
proximated by the sequential numerical integration of 
a one-particle system with complex boundary condi- 
tions. Similar as EMD, BTR is also much more ef- 
ficient than MD at the price of a restricted applica- 
bility, see Sec. [7J While the main precondition of 
EMD is vanishing contact time, for BTR we have to 
require that a particle having arrived in a stable po- 
sition does not leave this position anymore. That is, 
complementary to EMD, BTR requires persistant con- 
tacts. There are cases when the application of BTR 
leads to unphysical results, however, if the algorithm 
is applicable, it leads to a vast increase of numerical 
efficiency; at present we can simulate systems of many 
millions particles on a desktop computer. The range 
of its applicability and benchmarks of the algorithm 
are discussed at the end of this article. 



2 Bottom-to-top Reconstruction 

The idea of BTR goes back to Visscher and Bol- 
sterli [36] who suggested an algorithm that allows for 
the fast simulation of large static granular packings. 
The fundamental idea of their method is to consider 
the motion of the particles sequentially, unlike in MD 
simulations, where the coupled system of Newton's 
equations of motion is solved for all particles simul- 
taneously. 

The deposition of the particles of a granular pack- 
ing, e.g., a heap on the plane (x,y, 0), proceeds 
as follows: the first particle is inserted at position 
(x 1 1 mt , yf 1 *, 2? ut ). The coordinate 2? ut should be larger 
than the expected final height of the heap, the coordi- 
nates x 1 1 mt and y| nlt can be chosen at random (particles 
are scattered over a certain area) or can be fixed, e.g., 
(a;™ 1 *, y™ 1 *) = (0, 0), to simulate the build-up of a heap 
from a point source. The particle then falls until it 
touches the ground at (x™ 1 * , yf 11 * , R\). At this position 
the particle remains fixed. Then the second particle 
is inserted at position (x 2 mt , y 2 mt , z 2 mt ). It f ans down 
until it touches either the ground at (x 2 mt , y 2 mt , R2) or 
the first particle, whatever happens first. If it touches 
the ground it remains fixed there just like the first 
particle. If it, however, touches the first particle, it 
rolls down its surface in the downslope direction un- 
til it either touches the ground (if R% > Ri) where it 
remains fixed or (if R2 < Ri) until it loses contact at 
Z2 = R\ when both centers are at the same height. 
From there it falls to the ground where it remains 
fixed. The next particles are treated likewise. A par- 
ticle remains fixed if it either touches the ground or if 
it attains a local minimum where it is supported by 
two (in two dimensions) or three (in three dimensions) 
other already fixed particles. 

Thus, each particle moves according to a set of rules 
from one state to the next. In this sense the algorithm 
belongs to the class of event-driven algorithms. 
In Fig. [T] the algorithm is sketched for a two-dimen- 
sional system. The moving particle is drawn with a 
dotted line, fixed particles are drawn with a solid line. 
The algorithm by Visscher and Bolsterli [36] was later 
improved and generalized, e.g., [17-19] and also ap- 
plied to certain dynamical systems [4,5], see below. 
The BTR algorithm is not exact, that is, it does not 
solve the coupled set of Newton's equations for the 
dynamics of the many-particle systems. Instead, there 
are three major simplifying assumptions 

1. The motion of each particle i is computed with 
the wall and the other already-deposited parti- 
cles j = 0, — 1 considered as fixed obsta- 
cles. The positions of the particles j = . . . i — 1 
are not influenced by the motion of the new 
particle %. Thus, the trajectory of each parti- 
cle is computed while taking gravity as the only 
driving force into account. The other particles 
j = . . . i — 1 and the wall establish the (compli- 
cated) boundary conditions to this motion. This 
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Particle 1 is inserted and moves downward until it touches the 
ground. 

( 2 ') 
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Particle 2 is inserted and moves downward until it touches par- 
ticle 1. Then it rolls on the surface of particle 1 until it touches 
the ground. 

;' 3 'i 




Particle 3 is inserted and moves downward until it touches par- 
ticle 2. It then rolls to the right until it touches particle 1. This 
position is unstable (not a local minimum), thus, the particle 
continues to roll on particle 1 until it touches the ground. 



4 




Particle 4 first touches particle 2 and starts to roll to the right 
until it touches particle 3. This position is stable (local mini- 
mum), thus, particle 4 is fixed at this position. 



(5) 




Particle 6 first touches particle 4 and rolls on it until it loses con- 
tact when the centers of both particles are at the same height. 
Then the particle falls again until it touches particle 3. It con- 
tinues rolling until it touches particle 5 too where it is fixed 
(position is stable since it is a local minimum). 



Figure 1: Sketch of the BTR algorithm 

way, the system of Newton's equation for the N- 
body system is decoupled and N single-particle 
equations are solved instead. 

2. Collisions of particle i with the wall or with other 
particles are assumed to be perfectly inelastic. 

3. The time dependence of the particles' motion is 
disregarded. 

These simplifications increase the efficiency of the sim- 
ulation significantly. Jullien et al. [20] were able to 
simulate packings of up to 10 8 particles using this al- 
gorithm. The simplifications restrict the applicability 
of BTR to a rather small class of problems, see Sec. 
[71 If applicable, however, BTR is a very efficient sim- 
ulation tool. 
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3 Implementation of the BTR algorithm 

The BTR algorithm allows for a very efficient numeri- 
cal simulation. We briefly explain the implementation 
for the simulation of a two-dimensional heap which is 
the simplest application of BTR and restrict ourselves 
to the simulation of a two-dimensional system of poly- 
disperse particles. The extension to three dimensions 
is straightforward. A more detailed explanation can 
be found in [32] and the full source code of our pro- 
gram is available at [1]. 

To compute the deposition of the nth particle, we do 
not need the positions of all n — 1 previously deposited 
particles. Instead, it is sufficient to know the positions 
of particles that can come into contact with particle n, 
i.e., particles at the surface of the heap. This simplifi- 
cation is always justified for two-dimensional systems, 
whereas for three-dimensional systems it is only justi- 
fied if the ratio between the largest and smallest radii 
of the particles in the system does not exceed the Apol- 
lonian ratio TA = Rmax/Rmin = V% / (2 — a/3) ~ 6.46. 
Otherwise the smaller particles can fall into the gaps 
between the larger particles and thus, in principle, in- 
teract with all particles of the heap. In this case, 
we cannot exclude any pair interactions a priori. A 
natural choice of radii within (i? mm , i? m ax) is to as- 
sume that the total mass of particles from the interval 
(R, R + dR) is constant regardless of R. Such a distri- 
bution can be obtained from equi-distributed random 
numbers z £ [0, 1) using the transformation [32,33] 

D -Rmin^max /-.x 

1 ~~ ~r — tttr — Tp — \ • w 

We store the indices of the particles at the surface 
in the container surface of type map<double , int> 
where the key is the their x-coordinate. 
Consider the deposition of particle n. If n is not 
in contact with any other particle it moves down- 
ward until it touches the ground or another particle. 
Apart from the ground plane all particles on the sur- 
face of the heap whose centers belong to the interval 
Xn + Rn + Rm&x) are potential contact 
partners, see Fig. [2j The container map allows to very 
efficiently iterate through these contact candidates. 




Figure 2: The centers of the contact candidates of the falling 
particle n (crossed circle) belong to the interval (x n — R n — 
R max 5 %n + Rn + -Rmax) (black filled circles). The other surface 
particles are drawn gray, particles which are not part of the 
surface are drawn hollow 

For all contact candidates i that are located below 
the particle to be deposited, we check if and at which 



height yy both particles touch, disregarding for the 
moment possible interference from other particles, 
that is, we check whether 

yy = Vi + V(Rn + Ri) 2 - (x { - x n ) 2 , (2) 

has a real solution and determine the maximum of 
this expression for all i (see Fig. [3]). To determine the 
desired contact point also a possible contact with the 
ground must be taken into account. Particle n is now 



o 




Figure 3: For each of the candidates A-D (drawn in black in 
Fig. [2]| we determine the vertical coordinate yy of a possible 
contact. The highest contact is with candidate C, therefore, C 
is the contact partner we are looking for 

moved downward to contact its partner particle or the 
ground. In the latter case, its new position is stable 
by definition and its deposition is accomplished. If, 
however, the particle meets another particle first, its 
new position cannot be stable. Instead it continues 
its motion by rolling down the surface of the contact 
partner until it comes into contact with the wall or 
with another particle, or it continues to fall down, as 
sketched in Fig. [TJ 

Similar as described above, again all possible candi- 
dates for the additional contact partner are deter- 
mined, i.e., all particles that may be contacted by par- 
ticle n while maintaining contact to its present part- 
ner p. The new particle rolls to the left if x n < x p 
or to the right otherwise. Candidates for the con- 
tact are particles whose x-coordinate is in the interval 
x p . . .x p ± (Rp + 2R n + -R max ), see Fig. 01 Again the 
index range of the candidate partners can be obtained 
easily from the map surface. 




R[p]+2R[n]+Rmax R[p]-2R[n]-Rmax 



Figure 4: The centers of the candidates for the second partner 
particle are located within the circle of radius r = R p + 2R n + 
J?max around x p . Depending on the relative position of the par- 
ticles n and p the circle can be reduced to a half circle. The 
particles represented by filled circles are the candidates 

We investigate the possible contact of particle n rolling 
down the surface of particle p with each candidate 
i separately, disregarding all other particles for the 
moment, that is, we seek real solutions = (a^,y^) 
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of the system 

(3) 

| ffi f"i | — Rn ~i~ Ri i 

for contacts between particles n, p and n, i and of the 
system 

\r n Tp\ — R n -\- Rp 

(4) 

Vn=Rn. 

for a possible contact between particles n and p and of 
n with the floor. If these systems have real solutions 
the new position of particle n is determined by the so- 
lution r % n = (x l n ,y l n ) or r% = {x™;y l n ) with the largest 
maximum vertical component, provided, y l n or respec- 
tively y™ is larger than y p . In this case we found a 
new position of particle n. We now check for stability 
(i.e. if the particle is supported both from the right 
of his center and from the left). If it is stable the de- 
position of the particle is finished, otherwise particle i 
assumes the role of particle p and particle n continues 
to roll on its surface. If no candidate is found before 
particle n reaches y n = y p , it loses contact with p and 
again falls vertically downward, see last line in Fig. [TJ 
This procedure is repeated until a stable position of 
particle n is found. 

When the stable position of n is found, the map 
surface is updated. First, n is recorded as a mem- 
ber of the map since it became part of the surface. 
Second, after depositing n, other particles may be 
screened, such that no further particles may come into 
contact with them. Consequently, these particles are 
not members of the surface anymore and are removed 
from surface, see Fig. [5J 

M¥h 

Figure 5: Left: surface particles before the deposition of the 
new particle (drawn gray). Center: a new particle is deposited. 
Right: corrected list of surface particles 

The description of the algorithm and the sketch of 
the implementation is now complete. Figure [6] shows 
snapshots of a growing heap. The particles which are 
part of the surface are drawn filled. 
Since particles in the process of deposition can contact 
only the particles at the surface, unlike as in MD the 
computational complexity of depositing a particle n 
does not grow with the number of already deposited 
particles, 0(n) but it increases with 0(n 1/2 ) in 2d or 
with 0(n 2//3 ) in 3d which is the reason for the com- 
putational power of BTR. We discuss the efficiency of 
the algorithm below in Sec. [71 




Figure 6: Snapshots of the growing heap. The number of parti- 
cles is N = 100 and TV = 1400 

4 Example: Stratification in a Sand Heap 

When a heap of particles is created by sequentially 
depositing particles of different size, size segregation 
(stratification) is observed, which is caused by dif- 
ferent angles of repose for large and small particles 
[10,22,37,38]. The particles form stripes as shown in 
Fig. [JJ In three dimensional dunes and ripples, other 




Figure 7: Formation of stripes in a heap consisting of particles of 
different properties (experimental result). The figure was taken 
from [15] 

more complex stratification patterns are observed [2]. 
This effect has been studied and modeled extensively, 
e.g., in [8,9,12-14,16,21,23-27,27-30]. Similar struc- 
tures have been observed in sand overflown by wind 
or water [35]. 

Figure [8] shows a heap of N = 10 6 particles of two 
different radii. The effect of stratification is visible in 
the close-ups. The small particles are drawn filled. 

5 Sedimentation of Nano-Powders 

The BTR-method can be generalized to clusters of 
contacting particles. Here clusters are considered as 
perfecly rigid, that is, the particles belonging to a 
cluster do not change their relative positions. This 
idealized behavior may be adequate for cohesive nano- 
powders where the attractive surface forces exceed by 
far inertial forces, due to their enourmous surface area 
per unit mass. The attracting force between particles 
in contact provides a mechanism which fixes particles 
to a position where they where first deposited. Also 
due to the surface forces, nano-powders are frequently 
very pourous. In this section we describe briefly the 
application of BTR to coarsening of nano-powders in 
the process of repeated siphoning the material from 
one container into another. For a more detailed de- 
scription see [34]. 
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Figure 9: Typical cluster after a few disposition cycles. Dark 
circles represent surface particles. The surface is not contiguous, 
the gaps are, however, too small for outside particles to touch 
inner particles. 





Figure 8: A heap consisting of N = 10 6 particles of two different 
radii. Size segregation (stratification) can be seen in the close- 
ups 



Figure 10: A heap of 6 million particles. Left: Particles de- 
posited in clusters of 3. Right: After 18 redeposition cycles. 



For clusters, the first step of the BTR algorithm is 
handled as for simple spheres: We determine the po- 
sition of the falling cluster when one of its particles 
touches a particle of the sediment or the wall. 
Rolling of the cluster on the already deposited heap 
particles is more complicated. Apart from situations 
discussed in Sec. [3]one is confronted with a large num- 
ber of special cases whose discussion is outside the 
scope of this work [34]. 

For efficient simulation of both processes, falling and 
rolling, we need the accessible surface of the cluster, 
that is, those particles which due to geometry can 
come into contact with other particles or the floor. 
There are several efficient algorithms to determine this 
surface [31]. Fig. [9] shows an example cluster and 
its surface (filled circles). The holes in the surface 
are small enough to keep outside particles away from 
the inner particles. In our simulation of a coarsen- 
ing nano-powder of N = 6 • 10 6 cohesive particles, 
R G (0.9, 1.1), initially the particles (or small clusters 
of a few particles) are placed into a rectangular con- 
tainer of length L = 8000 limited by vertical walls. 
After the initial sedimentation due to standard BTR 
they are densely packed in the container (see Fig 1 101 
left). To obtain a flat surface of the heap one has to 
insert the particles at random positions distribuited 
uniformly over the length of the box. The material is 
now cut into square blocks of about 50 x 50 average 
particle diameters. The blocks are decomposed into 



clusters of mutually contacting particles. Now these 
clusters are sedimented into the container. Figure [11] 
shows the average height (h) of the material surface 
over the number of sedimentation cycles. 




Figure 11: Average height of the heap vs number of disposition 
cycles. The solid line is the fit to an exponential function, the 
diamonds are the measured values. 

This quantity was computed by cutting the box into 
narrow vertical strips and determining the height h of 
the highest particle in each strip. The height of the 
strip is now h + R where R is the radius of the highest 
particle. Averaging over all strips yields (h). The 
data points fit an exponential law of the form (h) = 
9450 — 5000 exp(— n/6) with n being the cycle number. 
The inital height does not fit into this scheme. This 
is due to the fact that the initial heap is build up 
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from very small structures (clusters of size 3 in the 
example) while each later sedimentation cycle is done 
with structures of a typical size 50. Fig. [TUl shows the 
heap after 18 sedimentation cycles. 

6 Dynamic Simulations 

BTR cannot be applied directly to the simulation of 
dynamic processes since due to the main principle of 
this algorithm, deposited particles cannot leave their 
positions anymore. A very restricted class of dynami- 
cal systems can be simulated, however, if we partition 
the dynamics into alternating steps of collective mo- 
tion and sequential deposition [19]. 
For the example of granular flow in a partially filled, 
slowly rotating cylinder, the partition of the (contin- 
uous) dynamics is [4,5] (see Fig. [T2"j) : 

1. initialization: Place the particles at random in- 
side of the container. 

2. collective motion: The positions of all particles 
follow the motion of the container for a small 
time step At. 

3. Preparation of the current time step: Increase 
the y-coordinate of all particles by a constant, 
e.g., i? cy i/10, with R cy \ being the cylinder radius. 

4. BTR: Apply BTR to the particles in sequence of 
increasing y-coordinate. 

5. Loop: Increase the system time by At and con- 
tinue with step [2j 




Figure 12: Dynamic BTR for the example of a slowly rotating 
cylinder, (a) The system at time t. (b) The cylinder is rotated 
(angle of rotation appears exaggerated), (c) All particles are 
lifted by fl cy i/10. (d, e) BTR, (J) Situation at time t + At when 
all particles are deposited 

Figure H3] shows a snapshot of a simulation of N = 10 6 
particles of different radii Ri G (0.1, 1) cm in a cylinder 
of radius 70 cm. The radii are chosen randomly in such 
a way that the the total mass of all particles from the 
interval (R, R + dR) is constant regardless of R. We 
notice that the small particles are concentrated close 
to the center of the cylinder. This effect, which is 
observed also experimentally, was found in simulations 
in [4,6]. 




Figure 13: Snapshot of a slowly rotating cylinder filled with 
N — 10 6 granular particles and close-ups. After one revolution 
consisting of 100 deposition steps, size segregation is clearly vis- 
ible. For better presentation in the top figure only 30% of the 
particles are drawn 

7 Benchmarks of the Implementation and 
Critical Analysis of the Model 

BTR-simulations perform in general much faster than 
regular MD. Although BTR was successfully applied 
to various large granular systems, this method is not 
universal, and there are cases where physically incor- 
rect behavior is observed [3]. 

A deposited particle does not move under the influ- 
ence of particles deposited later, even if the parti- 
cle suffers (in the realistic system) violent collisions. 
The motion of the particles is, thus, not governed 
by Newton's equation of motion. Instead, each sin- 
gle grain performs an overdamped motion in a com- 
plicated potential landscape comprising the already 
deposited grains. This is equivalent to disregarding 
inertial forces and moments. Even from a very ba- 
sic and intuitive concept of classical mechanics it is 
clear that this algorithm cannot describe the granular 
many-body problem in a general way. It should be 
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regarded as a compromise between computing power 
requirements and realistic description of physical re- 
ality. Figure demonstrates that BTR may lead to 
non-physical descriptions. 




Figure 14: left: The configuration of particles is physically pos- 
sible amd stable, but cannot be generated by BTR since none of 
the gray particles is in a stable position without the others, i.e., 
none of them could be deposited first, middle and right: BTR 
may generated this sequence, although it is not realistic since 
the configuration on the right-hand side is unstable 

Nevertheless, i/BTR is applicable, it leads to a great 
increase of performance as compared with MD. Fig- 
ure [TS] shows the CPU time to deposit a heap if N 
particles, using a desktop computer (Intel Pentium 4, 
3GHz). The program code is available at [1]. 




4 5 6 

10 10 10 

number of particles N 



Figure 15: CPU time to build a heap of N particles using BTR. 
The data points (diamonds) are in close agreement with a pow- 
erlaw t oc N 3 ^ 2 (solid line) as discussed in the text 

The algorithm may yield satisfactory results when 
treating systems where the boundary conditions 
change only very slowly. In such cases MD is fre- 
quently inefficient. In particular, for problems where 
the computation of the trajectories of individual par- 
ticles is less important, the algorithm can be seen as a 
good compromise between efficiency and precision of 
the result. 
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